DALS019-统计模型2-贝叶斯分布与层次分析

前言

这一部分是《Data Analysis for the life sciences》的第7章统计模型的第2小节，这一部分的主要内容涉及贝叶斯统计的一些原理。相应的R markdown文档可以参考作者的Github，另外，如果想补充一些贝叶斯的相关知识，可以参考这本书《统计学关我什么事：生活中的极简统计学》，这是一本有关贝叶斯统计的科普书，公式不多。

贝叶斯统计

高通量数据的一个明显特点就是，虽然我们最终会报道一些特定的基因，但是我们还会观察到许多相关的结果。例如，我们检测会数千个基因或者是代表了蛋白结合位点的数千个峰图，或者是一些CpGs的甲基化水平。但是，我们这里使用的多数统计推断方法都是独立地处理每个特征值，并且几乎忽略了来自其它特征的数据。在这一部分里，我们将会了解到，如何通过对特征值的联合建模来进行统计。这里方法中使用最为广泛的就是层次模型(hierachical models)，我们会在后面的贝叶斯统计中进行解释。

贝叶斯定理

先来看一个案例，如果我们有一种检测手段来检测囊性纤维化(cystic fibrosis)。假设这个检测手段的精确程度为99%，我们可以使用下面的公式来表示：

$$\mbox{Prob}(+ \mid D=1)=0.99, \mbox{Prob}(- \mid D=0)=0.99$$

其中，$+$表示阳性结果，$D$表示检测的结果，其中$1$表示得病，，$0$表示不得病。

现在我们随机选择一个人进行检测，结果如果是阳性，那么这个人患病的概率是多大？也就是说要计算$\mbox{Prb}(D=1|+)$的结果。囊性纤维化的发病率是每1/3900，也就是说，$\mbox{Prob}(D=1|+)=0.0025$，为了计算出这个人患病的概率，我们就会使用到贝叶斯定理，贝叶斯定理公式如下所示：

$$\mbox{Pr}(A \mid B) = \frac{\mbox{Pr}(B \mid A)\mbox{Pr}(A)}{\mbox{Pr}(B)}$$

这个公式就可以应用到我们的案例中，如下所示：

$$\begin{align*} \mbox{Prob}(D=1 \mid +) & = \frac{ P(+ \mid D=1) \cdot P(D=1)} {\mbox{Prob}(+)} \\ & = \frac{\mbox{Prob}(+ \mid D=1)\cdot P(D=1)} {\mbox{Prob}(+ \mid D=1) \cdot P(D=1) + \mbox{Prob}(+ \mid D=0) \mbox{Prob}( D=0)} \end{align*}$$

换成实际数字，则如下所示：

$$\frac{0.99 \cdot 0.00025}{0.99 \cdot 0.00025 + 0.01 \cdot (.99975)} = 0.02$$

也就是说，虽然这种检测手段有99%的精度，但是一个人的检测结果如果是阳性，那么这个人得病的概率只有0.02。这似乎有点反直觉。其原因就是，我们必须要考虑，当我们随机选择一个人时，这个人患上这种疾病时非常罕见的可能性。为了说明这一个，我们随便使用Monte Carlo模拟来计算一下。

模拟

下面的模拟旨在帮助你能够以可视化的形式来理解贝叶斯定理。我们首选从一个总体中随机选择1500人，其中患病的概率是5%，代码如下所示：

set.seed(3)
prev <- 1/20
##Later, we are arranging 1000 people in 80 rows and 20 columns
M <- 50 ; N <- 30
##do they have the disease?
d<-rbinom(N*M,1,p=prev)

现在进行一项检测，这个检测的准确率是90%，如下所示：

accuracy <- 0.9
test <- rep(NA,N*M)
##do controls test positive?
test[d==1] <- rbinom(sum(d==1), 1, p=accuracy)
##do cases test positive?
test[d==0] <- rbinom(sum(d==0), 1, p=1-accuracy)

由于没有患病的人数要远远超过患病的人数，即使存在着极低的假阳性率，那么在检测为阳性的结果中，不患病人的也要多于患病的人，如下所示：

cols <- c("grey","red")
people <- expand.grid(1:M,N:1)
allcols <- cols[d+1] ##Cases will be red
positivecols <- allcols
positivecols[test==0] <- NA ##remove non-positives
negativecols <- allcols
negativecols[test==1] <- NA ##remove non-positives
library(rafalib)
mypar()
layout(matrix(c(1,2,1,3),2,2),width=c(0.35,0.65))
###plot of all people
plot(people,col=allcols,pch=16,xaxt="n",yaxt="n",xlab="",ylab="",
     main=paste0("Population: ",round(mean(d)*100),"% are red"))
plot(people,col=positivecols,pch=16,xaxt="n",yaxt="n",xlab="",ylab="",
     main=paste("Tested Positive:",round(mean(d[test==1])*100),"% are red"))
plot(people,col=negativecols,pch=16,xaxt="n",yaxt="n",xlab="",ylab="",
     main=paste("Tested Negative:",round(mean(d[test==0])*100,1),"% are red"))

现在解释一下上面的图：

顶图：红点表示患者。每个人接受检测的话，有90%的可能性结果是准确的。即$\mbox{Pr}(D=1)$；

下左图：检测结果为阳性（无论结果对不对，这里面含有真阳性与假阳性）的人，即$\mbox{Pr}(D=1|+)$

下右图：检测结果为阴性的人，即$\mbox{Pr}(D=0|+)$

贝叶斯的实际运用

José Iglesias是一名职业棒球运行员，在2013年4月份，他开始了职业生成，他表表现得很好，成绩如下所示：

Month	At Bats	H	AVG
April	20	9	.450

平均击球率(battingaverage, AVG)统计是一种检测成功的方式。粗略地说，这个指标是告诉我们击球时的成功情况。AVG=0.450就表明，这个运动员在他击球的时间内（对应上面的At Bats）成功了45%，这是一个非常高的水平。为什么水平高呢，因为自1941年Ted Williams的AVG=0.400以来，还没有人能超过这个水平。为了说明层次模型的强大功能，我们将会在后面对Jose的数据进行预测。

在本书的前面部分到此为止，我们所到到的统计学技术可以被称为频率学派技术(frequentist techniques)，使用这种知识做出的结论就是能计算出一个置信区间。我们可以把击打(hitting)这个事件看作是一个二元结果，其成功的概论为$p$。因此，如果成功率是0.450的话，那么击球次数是20次时，标准误为：

$$\sqrt{\frac{.450 (1-.450)}{20}}=.111$$

也就是说，我们计算出的置信区间为.450-.222 to .450+.222或.228 to .672。

使用这种手段进行预测存在着两个问题。第一，实际用处不大。第二，成功率在0.450之间波动，也就是说这个人打破了Ted William的纪录。如果你自己关注棒球运动的话，打破Ted William纪录这种描述似乎是有问题，这是因为你隐含地使用了一个层次模型，这个模型会影响后面几年棒球的信息。在这里，我们对这种直觉进行量化。

首选，我们来研究一下前三个赛季里所有超过500次击打(at bats)的运动员的击球率的分布情况：

tmpfile <- tempfile()
tmpdir <- tempdir()
download.file("http://seanlahman.com/files/database/lahman-csv_2014-02-14.zip",tmpfile)
##this shows us files
filenames <- unzip(tmpfile,list=TRUE)
players <- read.csv(unzip(tmpfile,files="Batting.csv",exdir=tmpdir),as.is=TRUE)
unlink(tmpdir)
file.remove(tmpfile)

library(dplyr)
library(rafalib)
mypar(1,3)
for(y in 2010:2012){
  dat <- filter(players,yearID==y) %>% mutate(AVG=H/AB) %>% filter(AB>500)
  hist(dat$AVG*1000,xlab="AVG",freq=FALSE,main=y,xlim=c(200,360))
}

我们注意到，普通运行员的AVG为0.275，总体运动员的标准差为0.027。所以，我们看到，0.45这个数字与它们假设偏离非常大，因为这个数字离平均值的差距超过了6个标准差。Jose的这个数字是由于运气，还是说他确实是过去50年来最好的运动员？或者是这两者的结合？但是，他有多幸运，以及运动天赋有多高？如果我们确信他是由于运气出现的这个数字，我们应该把他送到相信他确实是0.45这个水平的球队里，并且有可能高估了他的潜力。

层次模型

层次模型为我们提供了如何观察到0.45的这个数字的数学描述。首先，我们会随机选择一个内在能力为$\theta$运动员，然后，我们看到成功概率为$\theta$的20个随机结果（这一段不太懂，原文为：The hierarchical model provides a mathematical description of how we came to see the observation of .450. First, we pick a player at random with an intrinsic ability summarized by, for example, $\theta$, then we see 20 random outcomes with success probability $\theta$.），其中：

$$\begin{align*} \theta &\sim N(\mu, \tau^2) \mbox{ describes randomness in picking a player}\\ Y \mid \theta &\sim N(\theta, \sigma^2) \mbox{ describes randomness in the performance of this particular player} \end{align*}$$

我们要注意2个层次（这就是为什么要称为层次分析的原因）：

1. 运动员与运动员之间的变异；
2. 击球时，运气因素导致的变异。

在贝叶斯框架中，第一个水平称为先验分布(prior distribution)，第二个水平称为采样分布(sampling distribution)。

现在我们使用贝叶斯模型来计算Jose的数据。假设我们想预测构成他真实击球平均水平$\theta$的内在能力的话。使用层次模型就按下面的方法表示：

$$\begin{align*} \theta &\sim N(.275, .027^2) \\ Y \mid \theta &\sim N(\theta, .111^2) \end{align*}$$

我们现在就可以计算出一个后验分布(posterior distribution)来描述我们对$\theta$的预测。这里我们可以使用贝叶斯规则的连续计算方法推导后验概率，后验概率是对给定观测数据参数$\theta$的分布：

$$\begin{align*} f_{ \theta \mid Y} (\theta\mid Y) &= \frac{f_{Y\mid \theta}(Y\mid \theta) f_{\theta}(\theta) }{f_Y(Y)}\\ &= \frac{f_{Y\mid \theta}(Y\mid \theta) f_{\theta}(\theta)} {\int_{\theta}f_{Y\mid \theta}(Y\mid \theta)f_{\theta}(\theta)} \end{align*}$$

我们主要是对能够使后验概率$f_{\theta\mid Y}(\theta\mid Y)$的值最大的$\theta$感兴趣。在我们的案例中，我们可以看出后验概率服从正态分布，我们能计算出均值$\mbox{E}(\theta\mid y)$，方差$\mbox{var}(\theta\mid y)$，尤其，我们可以计算出这个分布的均值服从以下分布：

$$\begin{align*} \mbox{E}(\theta\mid y) &= B \mu + (1-B) Y\\ &= \mu + (1-B)(Y-\mu)\\ B &= \frac{\sigma^2}{\sigma^2+\tau^2} \end{align*}$$

这是一个总体均值$\mu$和观测数据$Y$的加权均值。其权重取决于总体$\tau$的SD和我们观测数据$\sigma$的SD。这个加权均值有时候也会被称为shrinking，因为它缩小(shrink)了对先验均值的估计，在Joes Iglesias的数据中，结果如下所示：

$$\begin{align*} \mbox{E}(\theta \mid Y=.450) &= B \times .275 + (1 - B) \times .450 \\ &= .275 + (1 - B)(.450 - .275) \\ B &=\frac{.111^2}{.111^2 + .027^2} = 0.944\\ \mbox{E}(\theta \mid Y=450) &\approx .285 \end{align*}$$

方差如下所示：

$$\mbox{var}(\theta\mid y) = \frac{1}{1/\sigma^2+1/\tau^2} = \frac{1}{1/.111^2 + 1/.027^2} = 0.00069$$

标准差因此是0.026。我们开始时，使用了传统频率学派的思路计算出的95%置信区间忽略了来自于其他运动员的数据，因此会单纯地认为Joes的数据是0.0450 ± 0.220。我们随后使用了贝叶斯方法，整合了来源于其他运动员的数据，以及前几年的数据，计算出了后验概率。这种计算思路实际上就是经验贝叶斯方法，因此我们使用了数据先构建了先验知识。从后验结果中我们可以知道，通过报告一个以均值为中心的区间来报告所谓的95%的可信区间，其发生的概率为95%，在这个案例中，其结果是0.285 ± 0.052。原文：From the posterior we can report what is called a 95% credible interval by reporting a region, centered at the mean, with a 95% chance of occurring. In our case, this turns out to be: 0.285 ± 0.052.

贝叶斯可信区间表明，如果其他的球队发现了0.45这个数字，我们应该考虑到Joses可能转会到其他球队，因为我们预测到了Jose的水平高于平均不水平。有意思的是，Red Sox在7月份的时候将Jose转会到了Detroit Tigers队。这里是Jose在接下来的5个月内的击球率：

Month	At Bat	Hits	AVG
April	20	9	.450
May	26	11	.423
June	86	34	.395
July	83	17	.205
August	85	25	.294
September	50	10	.200
Total w/o April	330	97	.293

虽然这两个区间都包括了最终的击球平均值，但是贝叶斯可信区间提供了更精确的预测，尤其是，这种方法预测到Jose在本赛季的剩余时间里表现不佳。

练习

P308

层次模型

有关层次模型的内容可以参考作者的Github。

在这一部分里，我们会使用数据理论来描述在高通量数据分析中常用的方法。常规的思路就是构建一个两层的层次模型。一层用于描述样本/实验单元之间的变异，另外一层用于描述特征值之间的变异。这种分析方法类似于我们前面讲的棒球案例，即第一层用于描述不同运动员之间的变异，第二层用于描述一个运动员成功的随机性。我们这里𢪮的所有模型与方法都考虑了第一个变异水平，例如构建t检验的醋。第二个水平允许我们通过从所有的特征值里“借用(borrow)”信息用于对特征值进行统计推断，从而提供检验效能。

现在我们来看一个在基因表达数据中使用最为广泛的统计学方法。这个统计学方就是由limmaBioconductor包提供的。这个方法已经被用于改造分析RNAseq数据，例如edgeR《edgeR: a Bioconductor package for differential expression analysis of digital gene expression data.》和DESeq2《Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2》。这两个包提供了t检验的替代方案，它们通过对方差进行建模从而极大地改善了统计功效。然而在棒球案例中，我们是对均值进行了建模，这是与那两种方法建模的不同之处。对方差建模需要更深的数学知识，但是思路是一样的。我们以一个案例来说明一下这种方法。

下图是一个火山图，它显示了使用t检验来分析数据的结果，显示了效应大小(effect size)和p值，其中使用了6个重复样本（对照组3个，干预组3个），其中有16个基因是人为设定的差异基因。只有这16个基因的备选假设为真，在火山图上它们标记为蓝色，代码如下所示：

library(SpikeInSubset) ##Available from Bioconductor
data(rma95)
library(genefilter)
fac <- factor(rep(1:2,each=3))
tt <- rowttests(exprs(rma95),fac)
smallp <- with(tt, p.value < .01)
spike <- rownames(rma95) %in% colnames(pData(rma95))
cols <- ifelse(spike,"dodgerblue",ifelse(smallp,"red","black"))
with(tt, plot(-dm, -log10(p.value), cex=.8, pch=16,
xlim=c(-1,1), ylim=c(0,4.5),
xlab="difference in means",
col=cols))
abline(h=2,v=c(-.2,.2), lty=2)

上图是使用t检验计算了两组样本的差异基因，其中Spiked-in基因使用蓝色。剩下的基因中，p小于的用红色标明。

在上面的火山图中，我们将y轴的截止值(cut-off)设为了4.5，但是有一个蓝点的p值小于$10^{-6}$。但是，从这张图中我们会发现2点怪异之处。第一，按照5% FDR的标准，只有一个阳性结果是显著的，如下所示：

sum( p.adjust(tt$p.value,method = "BH")[spike] < 0.05)

结果如下所示：

> sum( p.adjust(tt$p.value,method = "BH")[spike] < 0.05)
[1] 1

这个结果与每组3个样本的低统计效能有关。第二，如果我们忽略掉统计推断，仅仅是基于t检验统计量的大小简单地对这些基因进行排序，那么我们会在任何大于1的排序列表中得到很多假阳性结果。例如，按照t检验统计量进行排序，位列前10名的基因中，有6个都是假阳性，如下所示：

table( top50=rank(tt$p.value)<= 10, spike) #t-stat and p-val rank is the same

结果如下所示：

> table( top50=rank(tt$p.value)<= 10, spike) #t-stat and p-val rank is the same
       spike
top50   FALSE  TRUE
  FALSE 12604    12
  TRUE      6     4

在火山图中，我们注意到，大多数基因的效能大小都非常小，这说明，估计的标准误非常小，我们可以通过画图的手段来看一下：

tt$s <- apply(exprs(rma95), 1, function(row)
  sqrt(.5 * (var(row[1:3]) + var(row[4:6]) ) ) )
with(tt, plot(s, -log10(p.value), cex=.8, pch=16,
              log="x",xlab="estimate of standard deviation",
              col=cols))

在这里我们就可以看到层次模型的用处了。如果我们假设这些变异的分布在所有基因中，然后我们通过分布来“调整”那些“太小”估计值，就可以改善我们的计算结果。在本书的前面部分中，我们提到F分布与观测到的方差分布近似，即：

$$s^2 \sim s_0^2 F_{d,d_0}$$

因为我们有数千个数据点，我们实际上可以检验一下这个假设，并且估计出参数$s_{0}$和$d_{0}$。这种估计方法指的就是经验贝叶斯统计，因为我们使用现有的数据（经验）就可以构建先验分布（贝叶斯方法）。

现在我们将前面的棒球案例应用到标准误的估计中。像以前一样，我们已经有了每个基因的观测值$s_{g}$，这是一个采样分布，用它来作为先验分布。我们因此可以计算出方差$\sigma_{g}^2$的后来又做分布，并且获得一个后验均值，细节可以参考文献《Linear models and empirical bayes methods for assessing differential expression in microarray experiments.》，均值如下所示：

$$\mathrm{E}\left[\sigma_{g}^{2} | s_{g}\right]=\frac{d_{0} s_{0}^{2}+d s_{g}^{2}}{d_{0}+d}$$

与棒球案例一样，后验均值会降低我们观测到的方差$s_{g}^2$偏向于全局方差$s_{0}^2$，其权重取决于样本大小，以及含有自由度$d$的样本数目，在这个案例中，就是取决于通过$d_{0}$的先验分布形状。（原文：AAs in the baseball example, the posterior mean _shrinks_ the observed variance $s_g^2$ towards the global variance $s_0^2$ and the weights depend on the sample size through the degrees of freedom $d$ and, in this case, the shape of the prior distribution through $d_0$. ）

在上面的图形中，我们可以看到40个基因的方差估计是如何缩小(shrink)的：

library(limma)
fit <- lmFit(rma95, model.matrix(~ fac))
ebfit <- ebayes(fit)
n <- 40
qs <- seq(from=0,to=.2,length=n)
idx <- sapply(seq_len(n),function(i) which(as.integer(cut(tt$s^2,qs)) == i)[1])
idx <- idx[!is.na(idx)]
par(mar=c(5,5,2,2))
plot(1,1,xlim=c(0,.21),ylim=c(0,1),type="n",
     xlab="variance estimates",ylab="",yaxt="n")
axis(2,at=c(.1,.9),c("before","after"),las=2)
segments((tt$s^2)[idx],rep(.1,n),
         ebfit$s2.post[idx],rep(.9,n))

上图显示的就是估计值如何向先验期望缩小的，40个基因包括了我们选择值的整个范围。

这种调整的一个重要方面就是使那些样本标准差接近于0的基因的样本偏差不再接近于0（向$s_{0}$收缩）。我们现在就创建一个t检验的统计模型，用于替代使用这个后验均值或“收缩“(shrunken)后的方差估计值。我们称这种t检验模型为适度t检验（moderated t-test）。当我们使用适应t检验后从火山图上就能明显地看到其改进之处：

library(limma)
fit <- lmFit(rma95, model.matrix(~ fac))
ebfit <- ebayes(fit)
limmares <- data.frame(dm=coef(fit)[,"fac2"], p.value=ebfit$p.value[,"fac2"])
with(limmares, plot(dm, -log10(p.value),cex=.8, pch=16,
                    col=cols,xlab="difference in means",
                    xlim=c(-1,1), ylim=c(0,5)))
abline(h=2,v=c(-.2,.2), lty=2)

这个火山图显示的就是使用适度t检验比较两组的差异基因结果。Spiked-in基因用蓝色进行了标注。剩下的基因中，p值小于的用红色标注。

现在我们来看一下排列前10的基因中假阳性的数目，这个数目就降为了2，如下所示：

table( top50=rank(limmares$p.value)<= 10, spike)

结果如下所示：

> table( top50=rank(limmares$p.value)<= 10, spike)
       spike
top50   FALSE  TRUE
  FALSE 12608     8
  TRUE      2     8

练习

P315